# Auxiliary functions for the ordered compositional data model

#================================================================
# TRANSFORMATION OF ORDERED COMP DATA TO R SPACE AND BACK
#================================================================

logit=function(x){log(x/(1+(1e-10)-x))}
a=function(y){(1+exp(y))/exp(y)}

v.y=function(v){
v.star=rep(NA,length(v)-1)
v=v[v>0]
v=logit(v[-1]/v[-length(v)])
v.star[1:length(v)]=v
return(v.star)
}

y.v=function(y){
y=y[!is.na(y)]
A=-a(y)*diag(length(y))
A=cbind(0,A)
A=rbind(A,1)
diag(A)=1
y.star=solve(A)[,length(y)+1]
return(y.star)
}


#================================================================
# Log-likelihood and sampling from low-truncated negative binomial
#================================================================


### sample from truncated negative binomial on (M,infty)
rtnegbin=function(N,mu,dispersion,l.bound){
sample=rep(NA,N)
mu=rep(mu,length=N)
l.bound=rep(l.bound,length=N)
n=1:N
ind=n
while(length(ind)!=0){ ##rejection sampling
sample[ind]=rnbinom(length(ind),mu=mu[ind],size=dispersion)
ind=which(sample<l.bound)
}
return(sample)
}


### negative binomial density on (l.bound, infty) where l.bound=M is number of seats in a district

dtlnegbin=function(x,mu,dispersion,l.bound){
N=length(x)
mu=rep(mu,length=N)
l.bound=rep(l.bound,length=N)
ldensity=ifelse(x>=l.bound,dnbinom(x,mu=mu, size=dispersion,log=TRUE)-
pnbinom(l.bound-.1,mu=mu,size=dispersion,lower.tail=FALSE,log=TRUE),log(1e-300))
return(ldensity)
}





#==============================================================
#---     Sweep operator for conditional MVN sampling
#==============================================================

cov.f=function(x,Sigma){
if (x[L+1]==L){return(x[-(L+1)])}
if (x[L+1]<L){
l=1:x[L+1]
t=(length(l)+1):L
S=Sigma[t,t]-Sigma[t,l]%*%solve(Sigma[l,l])%*%Sigma[l,t]
m=c(x[l],x[t]%*%chol(S))
return(m)
}
}

mvm=function(x,Sigma,g){
if (x[T+L+1]==L){return(rep(0,L))}
if (x[T+L+1]<L){
l=1:x[T+L+1]
t=(length(l)+1):L
m=x[1:T]%*%g
M=m[t]+Sigma[t,l]%*%solve(Sigma[l,l])%*%(x[(T+1):(T+length(l))]-m[l])
return(c(rep(0,length(l)),M))
}
}


#=================================================================
#     CENTER VARIABLES
#=================================================================

cntr=function(x){x-mean(x,na.rm=TRUE)}

#=================================================================
#     MAKE LAGS BASED ON FACTORS
#=================================================================

lag=function(x, steps=1, strat=NULL){
lp=function(x) {
R=sapply(1:steps,function(steps) c(rep(NA,steps),x[-((length(x)-steps+1):length(x))]))
R=matrix(R,ncol=steps)
colnames(R)=paste(paste(names(x),"lag",sep=""),1:steps,sep="-")
return(R)
}
if (length(strat)>0){
temp=tapply(x,factor(strat),lp)
return(do.call("c",temp))
}
return(lp(x))
}


#=================================================================
#      CREATE DESIGN MATRIX FOR FIXED OR RANDOM EFFECTS
#=================================================================

rem=function(x, fix.ef=TRUE){
x=factor(x)
list=lapply(levels(x),function(y){matrix(1*(x==y),length(x),1)})
M=do.call("cbind",list)
colnames(M)=levels(factor(x))
if (fix.ef==FALSE){return(M)}
return(M[,-1])
}


